	PROGRAM DCARBEB_REACT 
	 
c	This program uses  a y12m to invert a complex matrix set.

c	This program calculates the transport coefficients and properties
c	for swarms in arbitrarily configured electric and magnetic fields
c	under steady state conditions.  

c	Last Updated: 1/6/99 	

	INTEGER nsp,nsh,nn,ny,nep

	PARAMETER (nsp=100,nsh=7,nn=2*(nsp+1)*(nsh+2)**2,
     +	ny=2*(nsp+2)*(nsh+2)*(5*nsh+(nsp+2)*(nsh+3)),nep=500)
		! nsp is the max no. of Sonine polynomials
		! nsh is the max no. of spherical harmonics
		! nn is max no of rows in coefficient matrix
		! ny is the max no. of elements in coef matrix
	  
	INTEGER         ! Indices
     +  i,ij,j,jj,jjj,jjjmax,k,l,ld,lmax,lm,lx,
     +	m,md,mm,mmx,mdmx,mx,mb,mt,mmax,minimum,
     +  noe,nu,nud,numax,num,nux,iconv,nq

	INTEGER         ! General integer values
     +  dec,dist,kk,ndg,omr,n,mmm,nbt,spacing,inttemp,vnn

	DOUBLE PRECISION dfm,ll

	DOUBLE PRECISION alf,a(0:nsh,0:nsp,0:nsp),
     +  beta,bonn,bx,C1(0:nsp,0:nsp,0:nsh),
     +  C2(0:nsp,0:nsp,0:nsh),
     +  del(0:nsp,0:nsp,0:nsh,0:nsh,-nsh:nsh),        
     +  delta,e,eon2k,eonn,E0,J0,kb,me,
     +  monk,muon2k,Mn,n_o,p,pi,psi,psideg,sqrootoftwo,
     +  Nbar(0:nsp,0:nsp,0:nsh),
     +  rr,rrp,R0,R1,
     +  sigma0,Tb,To,temp,
     +  u,wi,wn,wmu,xi,zeta

	DOUBLE PRECISION         
c     +  gamlog(50),faclog(50),vx(15,25,25),as(25,25,6)
c     +  gamlog(90),faclog(90),vx(9,62,62),as(50,50,10)
     +  gamlog(150),faclog(150),vx(10,100,100),as(100,100,9)
c     +  gamlog(150),faclog(150),vx(9,82,82),as(82,82,9)
c     + gamlog(250),faclog(250),vx(10,200,200),as(200,200,9)	

!	DOUBLE PRECISION ! From collmatrix program
!     +  vrates(30,82),vrate(30,82),rate(30),
!     +  cwi,cw1,cw2,cTo,cTb
!	INTEGER ! From collmatrix program
!     +  cnumax,clmax

	DOUBLE PRECISION
c     + vrate0(2,-30:30,0:200),vrate1(2,-30:30,0:200),
     + vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100),
     + rate(2,-30:30),prate(2,-30:30),N_nu00,
     + field_power,
     +  cwi,cw1,cw2,cTo,cTb
	INTEGER ! From collmatrix program
     +  cnumax,clmax

	INTEGER
     + np(2)

	DOUBLE COMPLEX   A0,B0,
     +  BL(0:nsh,0:nsp,0:nsp,-nsh:nsh,-nsh:nsh),
     +  psi000(0:nsp,0:nsh,-nsh:nsh),
     +  psi11n1(0:nsp,0:nsh,-nsh:nsh),
     +  psi110(0:nsp,0:nsh,-nsh:nsh),psi11p1(0:nsp,0:nsh,-nsh:nsh),
     +  temporary,
     +  tmp,tmpa,tmpb,tmpc,tmp1,tmp2,tmp3,tmp4

	DOUBLE PRECISION ! Now when IMSL not used
     +  rhs000(nn),AA0(ny),
     +  rhs110(nn),rhs11p1(nn)

	DOUBLE PRECISION  ! Transport properties 
     +	en0,DXX,DXY,DXZ,DYX,DYY,DYZ,DZX,DZY,
     +  DZZ,gamx,gamy,gamz,gam,phi,theta,Te,
     +  Txx,Txy,Txz,Tyy,Tyz,Tzz,W_x,W_y,W_z,Wtot,
     +  omegaonn_o,W_xb,W_yb,W_zb,Wtotb

	DOUBLE PRECISION ! tmp quantities
     +  tmpa1,tmpa2,tmpa3,tmpa4,tmpa5,tmpa6,tmpa7,tmpa8,tmpa9,
     +  tmpdiffa,tmpdiffusion,
     +  tmpenergya,tmpenergyb,tmpenergy,
     +  tmpdrifta,tmpdriftb,tmpdrift,tmptempa	


c	FOR THE IMSL COMPLEX INVERSION ROUTINES

c	INTEGER         ! For inversion process
c    +  iparam0(6),ipath0,ipvt0(nn),irfac0(3*ny),nz0,irow0(ny),
c     +  jcol0(ny),jcfac0(3*ny),jpvt0(nn),nfac0,nl0

c	DOUBLE PRECISION            ! For inversion process
c     +  rparam0(5)

c	DOUBLE COMPLEX      ! For complex inversion process imsl
c     +  fac0(3*ny)

c	FOR THE Y12M REAL INVERSION ROUTINES

    	INTEGER ! For y12m inversion  
     +  nz0,nz1,nz2,ha0(nn,11),iflag0(10),iflag1(10),ifail0,
     +  irow0(ny),jcol0(ny)

	DOUBLE PRECISION 
     +  pivot0(nn),pivot1(nn),pivot2(nn),aflag0(8),aflag1(8),aflag2(8)

	DOUBLE PRECISION  ! For distribution function calculation
     +  sonine,
     +  en,c,cc,delen(0:6),sv,svv,energy(nep),x,
     +  enx,eny,costheta,ALP(0:nsh,0:nsh),phiplane,vdfspacing,
     +  distfn(0:90,0:nep,0:nep),disti,
     +  xen(0:nep),yen(0:nep),xx,delenx,deleny

	DOUBLE COMPLEX
     +	f0i,f0(0:nsh,0:nsh,nep)

	INTEGER noexp,noeyp,
     +  enptsi,enpts(0:6),pp,dist1

c	COMMON/WEALTH/FACLOG,GAMLOG
	COMMON/WEALTH/FACLOG,GAMLOG,VX,AS
c	EXTERNAL dlftzg,dlfszg,dlslzg

	CHARACTER filename*30,datetime*24,lable*40

c        DATA delen/0.000001d0,0.001d0,0.01d0,0.02d0,0.05d0,0.1d0,0.5d0/
c        DATA enpts/1,10,9,15,12,5,11/
  	DATA delen/0.05d0,0.02d0,0.02d0,0.05d0,0.1d0,0.2d0,0.5d0/
        DATA enpts/500,50,50,45,10,5,10/

c*****  Defining Input Parameters  **********************************

	open(unit=4,file='dctwotemp.dat',status='old')
	open(unit=2,file='eonn.dat',status='old')
	open(unit=3,file='Tb.dat',status='old')
	print*, 'ny=', ny
        print*, 'Input the electric field amplitude (Td)'
	read(2,*) eonn
	print*, eonn
       print*, 'Input the magnetic field amplitude (Hx)'
	read(4,*) bonn
	print*, bonn
       print*, 'Input the angle between E and B (psi)'
	read(4,*) psideg
	print*, psideg
       print*, 'Input the ion mass (amu)'
	read(4,*) wi
	print*, wi
       print*, 'Input the neutral mass (amu)'
	read(4,*) wn
	print*, wn
       print*, 'Input the background temperature T_o (K)'
	read(4,*) To
	print*, To
       print*, 'Enter the max,min number of Sonine polynomials desired'
	read(4,*) nux,num
	print*, nux,num
       print*, 'Enter the max,min number of spherical harmonics desired'
	read(4,*) lx,lm
	print*, lx,lm
	print*, 'Enter the max,min number of m indices'
	read(4,*) mx,mm
	print*, mx,mm	
       print*, 'Enter the order of the density gradient expansion'
	read(4,*) ndg
	print*, ndg
       print*, 'Enter the order of the mass ratio expansion'
	read(4,*) omr
	print*, omr
       print*, 'Do you the distribution function outputed, 0=yes'
	read(4,*) dist
	print*, 'Enter the plane of the dist to be calculated'
	read(4,*) phiplane
	print*, 'Enter the spacing between vdf calculations'
	read(4,*) vdfspacing
	print*,dist
       print*, 'Power Law Cross section of the form bc^p'
       print*, 'Enter b and p'
	read(4,*) bx,p
	print*, bx,p
	print*, 'Enter the basis temperature'
	read(3,*) Tb
	print*, 'Tb=',Tb
c       print*, 'Model Used Label'
	read(4,*) lable
c	print*, 'Enter the degree spacing'
	read(4,*) spacing

c*****  Defining constant parameters used ************************************

c       Source: Physics Today, August 1990, Part2, Seventh Annual Buyers Guide

	kb= 1.380658d-23             !Boltzmann's constant (JK^-1)
	e=1.60217733d-19           !electronic charge (C)
	u=1.6605402d-27             ! 1 amu (kg)
	me=5.48579903d-4            ! Electron mass in amu
	pi=3.141592653589793238d0
	pi=4.d0*datan(1.d0)
	psi=2.d0*pi*psideg/360.d0
	

c       (This routine reads in data and calculates faclog,gamlog and
c        massratioinputs if required)

        print*, 'into data'
	print*, '*****=',lx,nux
c        call datainput(To,wi,wn,lx,nux)
	print*, '*****=',lx,nux
        print*, 'outof data'
c	print*, 'wn=',wn
c	wn=4.d0
c	wn=16.d0

	wmu = (wi*wn)/(wi+wn)*u   !Reduced mass
	eon2k = e/2.d0/kb         ! e/2k
	monk= wi/kb*u           ! m/2k
	muon2k = wmu/2.d0/kb    ! mm_o/(m+m_o)2k
	Mn=wn/(wn+wi)           ! Something like reduced mass of neutrals
	sigma0 = 1.d-20         ! Scaling cross section 1A^2
	sqrootoftwo=dsqrt(2.d0)

	alf=dsqrt(monk/Tb)            ! alpha parameter 
	zeta=0.1d0*dsqrt(Mn)*eon2k/Tb*eonn
					! Dimensionless electric field term
	beta=1d-7*dsqrt(muon2k/Tb)*e/me/u*bonn
					! Dimensionless magnetic field term
	print*, 'beta=', beta
	print*, 'zeta=', zeta

c       (Making array faclog (log of factorial of intergers)
c       (faclog(i)=dlog((i-1)!)

	faclog(1)=0.d0
	do 919 i=2,nux+4
		xi=dfloat(i-1)
		faclog(i)=faclog(i-1)+dlog(xi)
919     continue
	
c       (Making array gamlog (log of gamma function of half integers))
c       (Gamlog(i)=dlog(Gamma((2*i-1)/2)) )

	gamlog(1)=0.572364942924700176d0  ! = ln(dsqrt(pi))
	do 921 i=2,nux+lx+6
		xi=dfloat(2*i-3)/2.d0
		gamlog(i)=gamlog(i-1)+dlog(xi)
921     continue

	call fdate(datetime)
!	datetime='Greifswald 2002'
	write(36,222)datetime,lable
	write(35,222)datetime,lable
	write(34,222)datetime,lable
	write(33,222)datetime,lable
	write(32,222)datetime,lable
	write(31,222)datetime,lable
	write(30,222)datetime,lable
	write(24,222)datetime,lable
	write(23,222)datetime,lable
	write(22,222)datetime,lable
	write(21,222)datetime,lable
222     format(1x,20a,5x,80a)
	write(30,*)
	write(31,*)
	write(32,*)
	write(33,*)
	write(34,*)
	write(35,*)
	write(36,*)
	write(21,*)
	write(22,*)
	write(23,*)
	write(24,*)
	write(30,223)eonn,bonn,To,wn
	write(31,223)eonn,bonn,To,wn
	write(32,223)eonn,bonn,To,wn
	write(33,223)eonn,bonn,To,wn
	write(34,223)eonn,bonn,To,wn
	write(35,223)eonn,bonn,To,wn
	write(36,223)eonn,bonn,To,wn
	write(21,223)eonn,bonn,To,wn
	write(22,223)eonn,bonn,To,wn
	write(23,223)eonn,bonn,To,wn
	write(24,223)eonn,bonn,To,wn
	write(30,533)Tb
	write(31,533)Tb
	write(32,533)Tb
	write(33,533)Tb
	write(34,533)Tb
	write(35,533)Tb
	write(35,533)Tb
	write(36,533)Tb
	write(21,533)Tb
	write(22,533)Tb
	write(23,533)Tb
	write(24,533)Tb


223     format(1X,'E/n_o =',F9.6,' Td.',2X,'B/n_o=',f7.1,1x,'Hx',
     +  2x,'T_o = ',f5.1,1X,'K',2x,'M = ',1pd10.4,' amu')           
533	format(1x,'Tb= ',f9.1,2x,'K')
	write(30,*)
	write(31,*)
	write(32,*)
	write(33,*)
	write(34,*)
	write(35,*)
	write(36,*)
	write(30,225)
	write(32,225)
225     format(1x,'l',1x,'m',1x,'nu',5x,'E',11X,'WX',11X,'WY',11X,
     +   'WZ',11x,'W',3x,'R/n_o')
	write(31,226)
226     format(1x,'l',1x,'m',1x,'nu',5x,'DXX',9X,'DYY',9X,'DZZ',8x,
     $            'DXZ+DZX',5X,'DXY+DYX',4X,'DYZ+DZY')
	write(33,327)
327     format(1x,'l',1x,'m',1x,'nu',5x,'DXZ',9x,'DZX',9x,'DXY',9x,'DYX',
     $	9x,'DYZ',9x,'DZY')
c	write(33,328)
c328     format(21x,'DYZ',12x,'DZY')
	write(34,329)
329     format(1x,'l',1x,'m',1x,'nu',3x,'gamx',8X,'gamy',8X,'gamz',8x,
     $             'gam',7x,'theta',2x,'phi')
	write(35,330)
330     format(1x,'l',1x,'m',1x,'nu',5x,'Txx',9x,'Tyy',9x,'Tzz',9x,'Txy',
     $	9x,'Txz',9x,'Tyz',9x,'Te')
	write(36,331)
331     format(1x,'l',1x,'nu',4x,'E',8x,'W',8x,'DL',8x,'DT',
     $	8x,'T_L',8x,'T_T',8x,'Gam')
	write(30,*)
	write(31,*)
	write(32,*)
	write(33,*)
	write(34,*)
	write(35,*)
	write(36,*)
	write(21,*)
	write(22,*)
	write(23,*)
	write(24,*)

c------ Evaluating the electric field reduced matrix elements ----------------

	 do 107 nu=0,nux		
	  do 107 nud=0,nux
	   do 107 l=0,lx
	    do 107 ld=0,lx
	     do 107 m=-l,l 
	    tmp=alf*dsqrt(2.d0/dfloat(2*l+1))*
     $    (dsqrt(dfloat(l-m)*dfloat(l+m)*(dfloat(nu+l)+0.5d0)/
     $    dfloat(2*l-1))*delta(nud,nu)*delta(ld,l-1) - 
     $    dsqrt(dfloat(l-m+1)*dfloat(l+m+1)*dfloat(nu)/
     $    dfloat(2*l+3))*delta(nud,nu-1)*delta(ld,l+1)) 

	   del(nu,nud,l,ld,m)=tmp
107	 continue
	      
c------ Evaluating the collision matrix at all Tbi --------------------------

c	     call coll(lx,nux,a,To,Tb)
c          call coll(Tb,To,Wn,Wi,lx+1,nux+1,p,bx)
c          call coll(Tb,To,Wn,Wi,lx+1,nux+1)

c	Checking the output from collmat and the parameters used
c	in this part of the code coincide.

	read(13) clmax,cnumax,cTo,cTb,cwi,cw1,cw2

	print*, 'Max lmax from coll = ',clmax-1,lx
	if (lx.gt.clmax-1) then
	    print*, 'lmax is inconsistent with collmat'
	    pause
	    stop
	endif
        print*, 'Max numax from coll = ',cnumax-1,nux
	if (nux.gt.cnumax-1) then
	    print*,'numax is inconsistent with collmat'
	    pause
	    stop
	endif
        print*, 'To from coll = ',cTo,To
	if (cTo.ne.To) then
	    print*,'To is inconsistent with collmat'
	    pause
	    stop
	endif
        print*, 'Tb from coll = ',cTb,Tb
	if (cTb.ne.Tb) then
	    print*,'Tb is inconsistent with collmat'
	    pause
	    stop
	endif
        print*, 'Ion mass from coll = ',cwi,wi
        print*, 'Neutral mass1 from coll = ',cw1,wn
	if (cw1.ne.wn) then
	    print*,'Masses are inconsistent with collmat'
	    stop
	endif
        print*, 'Neutral mass2 from coll = ',cw2

	read(14) as

		do 110 l=0,lx
		 do 110 nu=0,nux
                  do 110 nud=0,nux
	    a(l,nu,nud)=as(nu+1,nud+1,l+1)  
110	        continue

	read(12) np(1),np(2),vrate0,vrate1
        

         do 1090 nu=0,nux	
          do 1090 nud=0,nux
           do 1090 l=0,lx
             tmp=dsqrt(dexp(faclog(nu+1))*
     $          dexp(gamlog(nud+l+2))/
     $          dexp(faclog(nud+1))/dexp(gamlog(nu+l+2)))
             Nbar(nu,nud,l)=tmp
1090    continue

c------ Evaluating the Velocity reduced matrix element --------------------

	   do 116 nud=0,nux
	    do 116 nu=0,nux
	     do 116 l=0,lx
	if (l.eq.0) then
	  C1(nu,nud,0)=0.d0
	else
	  C1(nu,nud,l)=dsqrt(2.d0/dfloat(2*l+1))*
     $ 		(dsqrt(dfloat(nu+l)+0.5d0)*delta(nu,nud)
     $		-dsqrt(dfloat(nu+1))*delta(nud,nu+1))
	endif
116	continue

	   do 117 nud=0,nux
	    do 117 nu=0,nux
	     do 117 l=0,lx
	if (l.eq.lx) then
	 C2(nu,nud,l)=0.d0
	elseif (nu.eq.0) then
	C2(nu,nud,l)=dsqrt(2.d0/dfloat(2*l+1))*
     $		dsqrt(dfloat(nu+l)+1.5d0)*delta(nu,nud)
	else
	C2(nu,nud,l)=dsqrt(2.d0/dfloat(2*l+1))*
     $		(dsqrt(dfloat(nu+l)+1.5d0)*delta(nu,nud)-
     $ 		dsqrt(dfloat(nu))*delta(nud,nu-1))
	endif 
117	continue

c#####  Loop over all angles 0->pi/2 ###########################
	
	vnn=0	! Counter for vdf
c	phiplane=phiplane*pi/180.d0	! Convert to radians
c	temp=90.d0/dfloat(spacing)
c	jjjmax=int(temp)	
c	do 9999 jjj=0,jjjmax

c	psideg=jjj*dfloat(spacing)
	print*, 'psi= ',psideg
	psi=2.d0*pi*psideg/360.d0

	write(30,*) 'Angle= ',psideg
	write(31,*) 'Angle= ',psideg
	write(32,*) 'Angle= ',psideg
	write(33,*) 'Angle= ',psideg
	write(34,*) 'Angle= ',psideg
	write(35,*) 'Angle= ',psideg
	write(36,*) 'Angle= ',psideg
	write(21,*) 'Angle= ',psideg
	write(22,*) 'Angle= ',psideg
	write(23,*) 'Angle= ',psideg
	write(24,*) 'Angle= ',psideg
	write(30,*)
	write(31,*)
	write(32,*)
	write(33,*)
	write(34,*)
	write(35,*)
	write(36,*)
	write(21,*)
	write(22,*)
	write(23,*)
	write(24,*)

	do 178 nu=0,nux
	do 178 nud=0,nux
	 minimum=min(nu,nud)
	do 178  l=0,lx
	  do 178 m=-l,l
	    do 178 md=-l,l
	tmp=delta(nu,nud)*
     $  beta*((dsqrt(dfloat(l-m)*dfloat(l+m+1))*delta(md,m+1) 
     $ -dsqrt(dfloat(l+m)*dfloat(l-m+1))*delta(md,m-1))
     $  *0.5d0*dsin(psi)+
     $  dcmplx(0.d0,-1.d0)*dfloat(m)*dcos(psi)*delta(m,md))
	    BL(l,nu,nud,m,md)=tmp
178     continue

c***** Convergence Loop for l and nu ***********************************

c       To check the convergence in the l and nu indices we loop over
c       the range of lmax and numax. 

	do 999 lmax=lm,lx
	if (mx.gt.lmax) then
	  mt=lmax
	else
	  mt=mx
	endif
	if (mm.gt.mt) then 
      	  mb=mt
	else
	  mb=mm
	endif
	  do 999 mmax=mb,mt
	   do 998 numax=num,nux

	 print*, 'nu= ',numax
c------ Initializing the matrices ------------------------------------------

	 i=1
	 do 25 l=0,lx
	  do 24 m=-lx,lx 
	   do 20 nu=0,nux
		psi000(nu,l,m)=0.d0
		psi110(nu,l,m)=0.d0
		psi11p1(nu,l,m)=0.d0
		i=i+1
20         continue
24        continue
25       continue

c-----  Solving the (s,lambda)=(0,0) system of equations -------------------

	do 358 ij=1,2 		! Blocks for total matrix
	  i=0
	  do 358 l=0,lmax
	    if (l.ge.mmax) then
	     	mmx=mmax
	    else
	 	mmx=l
	    endif	
	    do 358 m=-mmx,mmx
	     do 358 nu=0,numax
		if ((l.eq.0).and.(nu.eq.0)) goto 358
		i=i+1
		rhs000(i)=0.d0
		rhs110(i)=0.d0
		rhs11p1(i)=0.d0
358	continue  
 	noe=i   ! No. of Complex equations


c......	Initial estimate of the reaction rate	....................

	if ((numax.eq.num).and.(lmax.eq.lm)) then
	 rrp=-Nbar(0,0,0)*a(0,0,0)
	endif
	iconv=0 	! convergence of eigenvalue

c ..... Building the matrices (s,lambda,mu)=(0,0,0) equations ............

1110	k=0		! Counter  for non-zero elements
	iconv=iconv+1
c	print*, 'iconv= ',iconv
	if (iconv.gt.50) then
	 print*, 'Too slow - convergence'
c	 stop
	 goto 1120
	endif

	i=1
	 do 35 l=0,lmax
	    if (l.ge.mmax) then
	     	mmx=mmax
	    else
	 	mmx=l
	    endif	
	  do 32 m=-mmx,mmx
	   do 30 nu=0,numax
	   j=1
	   do 29 ld=0,lmax
	    if (ld.ge.mmax) then
	     	mdmx=mmax
	    else
	 	mdmx=ld
	    endif
	    do 28 md=-mdmx,mdmx
	     do 27 nud=0,numax

	    if ((nu.eq.0).and.(l.eq.0)) goto 30 
	    if ((nud.eq.0).and.(ld.eq.0)) goto 27
				! i.e. wipe out first row
				! and place first column on rhs

	   E0=0.d0
	   if (m.eq.md) 
     $	     E0=-sqrootoftwo*zeta*del(nu,nud,l,ld,m)/alf

	   B0=dcmplx(0.d0,0.d0)
	  if ((l.eq.ld).and.(nu.eq.nud)) 
     $      B0=BL(l,nu,nud,m,md)	     

	   J0=0.d0
	   R0=0.d0
	   if ((l.eq.ld).and.(m.eq.md))  then
     	     J0=Nbar(nu,nud,l)*a(l,nu,nud)
	     if (nu.eq.nud) then
	       R0 = rrp
	     endif
	   endif

		A0=J0+E0+B0+R0

	     if (dabs(real(A0)).ge.(1.d-100)) then   ! Into sparse form
	      k=k+1
	      AA0(k)=real(A0)
	      irow0(k)=i
	      jcol0(k)=j
	      k=k+1
	      AA0(k)=real(A0)
	      irow0(k)=i+noe
	      jcol0(k)=j+noe
	     endif
	     if (dabs(dimag(A0)).ge.(1.d-100)) then
	      k=k+1
	      AA0(k)=-dimag(A0)
	      irow0(k)=i
	      jcol0(k)=j+noe
	      k=k+1
	      AA0(k)=dimag(A0)
	      irow0(k)=i+noe
	      jcol0(k)=j
	     endif
	    j=j+1
27          continue
28          continue
29         continue
	   i=i+1
30         continue
32        continue
35       continue
	nz0=k                           ! No. of non-zero elements in AA0
	kk=nz0

	i=1
	 do 50 l=0,lmax
	 if (l.ge.mmax) then
	  	mmx=mmax
	 else
	 	mmx=l
	 endif	
	 do 47 m=-mmx,mmx 
	  do 45 nu=0,numax
	   
           if ((l.eq.0).and.(nu.eq.0)) goto 45 

	   E0=0.d0
	   if (m.eq.0) 
     $	     E0=-sqrootoftwo*zeta*del(nu,0,l,0,m)/alf
	   J0=0.d0
	   if ((l.eq.0).and.(m.eq.0))  
     $      J0=Nbar(nu,0,0)*a(0,nu,0)

	    rhs000(i)=-real(J0+E0)
	    rhs000(i+noe)=0.d0 ! NB B0=0.d0 
	    i=i+1
45        continue
47       continue
50       continue

	mmm=2*noe  ! No. of equations for real inversion

	call initial(aflag0,aflag1,iflag0,iflag1)
	call y12mbf(mmm,nz0,AA0,jcol0,ny,irow0,ny,ha0,nn,aflag0,
     +               iflag0,ifail0)
	if (ifail0.ne.0) go to 5151
	call y12mcf(mmm,nz0,AA0,jcol0,ny,irow0,ny,pivot0,rhs000,
     +       ha0,nn,aflag0,iflag0,ifail0)
	if (ifail0.ne.0) go to 5151
 	call y12mdf(mmm,AA0,ny,rhs000,pivot0,jcol0,ha0,nn,iflag0,ifail0)
5151 	print*, 'Problem here in y12m',ifail0
	continue

	i=1
	do 52 l=0,lmax
	 if (l.ge.mmax) then
	  	mmx=mmax
	 else
	 	mmx=l
	 endif	
	 do 52 m=-mmx,mmx
	   do 51 nu=0,numax	
	    if ((l.eq.0).and.(nu.eq.0)) goto 51
	      psi000(nu,l,m)=rhs000(i)+dcmplx(0.d0,1.d0)*rhs000(i+noe) 
				     ! Converting back to standard form
	      i=i+1
51         continue
52      continue
	psi000(0,0,0)=1.d0

	rr=0.d0   ! New estimate of the reaction rate 
	do 531 nud=0,numax
	  rr=rr-Nbar(0,nud,0)*a(0,0,nud)*psi000(nud,0,0)
531	continue

	if (dabs(1.d0-(rr/rrp)).gt.(1.d-5)) then
	 rrp=rr
	 goto 1110
	else
	 goto 1120
	endif
1120	continue

c ----- Solving the first order density gradient equations ----------------

	if (ndg.ge.1) then

c	do 350 i=1,nz0  ! Initilaizing those element used in spatially uniform case
c	 AA0(i)=0.d0
c	 irow0(i)=0
c 	 jcol0(i)=0
c350	continue

	i=1
	k=0
	 do 351 l=0,lmax
	    if (l.ge.mmax) then
	     	mmx=mmax
	    else
	 	mmx=l
	    endif	
	  do 321 m=-mmx,mmx
	   do 301 nu=0,numax
	   j=1
	   do 291 ld=0,lmax
	    if (ld.ge.mmax) then
	     	mdmx=mmax
	    else
	 	mdmx=ld
	    endif
	    do 281 md=-mdmx,mdmx
	     do 271 nud=0,numax

	    if ((nu.eq.0).and.(l.eq.0)) goto 301 
	    if ((nud.eq.0).and.(ld.eq.0)) goto 271
				! i.e. wipe out first row
				! and place first column on rhs

	   E0=0.d0
	   if (m.eq.md) 
     $	     E0=-sqrootoftwo*zeta*del(nu,nud,l,ld,m)/alf

	   B0=dcmplx(0.d0,0.d0)
	  if ((l.eq.ld).and.(nu.eq.nud)) 
     $      B0=BL(l,nu,nud,m,md)	     

	   J0=0.d0
	   R0=0.d0
	   if ((l.eq.ld).and.(m.eq.md))  then
     	     J0=Nbar(nu,nud,l)*a(l,nu,nud)
	     if (nu.eq.nud) then
	       R0 = rr		! Using the solution from (s,lambda,mu)=(0,0,0)
	     endif
	   endif

	   R1=0.d0
	   if ((ld.eq.0).and.(md.eq.0))
     $	    R1=-Nbar(0,nud,0)*a(0,0,nud)*psi000(nu,l,m)

		A0=J0+E0+B0+R0+R1

	     if (dabs(real(A0)).ge.(1.d-100)) then   ! Into sparse form
	      k=k+1
	      AA0(k)=real(A0)
	      irow0(k)=i
	      jcol0(k)=j
	      k=k+1
	      AA0(k)=real(A0)
	      irow0(k)=i+noe
	      jcol0(k)=j+noe
	     endif
	     if (dabs(dimag(A0)).ge.(1.d-100)) then
	      k=k+1
	      AA0(k)=-dimag(A0)
	      irow0(k)=i
	      jcol0(k)=j+noe
	      k=k+1
	      AA0(k)=dimag(A0)
	      irow0(k)=i+noe
	      jcol0(k)=j
	     endif
	    j=j+1
271          continue
281          continue
291         continue
	   i=i+1
301         continue
321        continue
351       continue
	nz0=k                           ! No. of non-zero elements in AA0

	i=1
	 do 500 l=0,lmax
	 if (l.ge.mmax) then
	  	mmx=mmax
	 else
	 	mmx=l
	 endif	
	 do 470 m=-mmx,mmx 
	  do 450 nu=0,numax

	    if ((l.eq.0).and.(nu.eq.0)) goto 450 

	  tmp1=dcmplx(0.d0,0.d0)
	  tmp2=dcmplx(0.d0,0.d0)
	  do 442 nud=0,numax

	if (l.eq.0) then 
	tmp1=tmp1+((
     $	C1(nu,nud,l)*dsqrt(dfloat(l+m)*dfloat(l-m)/
     $	dfloat(2*l-1))*0.d0+
     $	C2(nu,nud,l)*dsqrt(dfloat(l-m+1)*dfloat(l+m+1)/
     $	dfloat(2*l+3))*psi000(nud,l+1,m))/alf)
	
	tmp2=tmp2+((
     $	C1(nu,nud,l)*dsqrt(dfloat(l+m)*dfloat(l+m-1)/2.d0/
     $	dfloat(2*l-1))*0.d0-
     $	C2(nu,nud,l)*dsqrt(dfloat(l-m+1)*dfloat(l-m+2)/2.d0/
     $	dfloat(2*l+3))*psi000(nud,l+1,m-1))/alf)
	else
	tmp1=tmp1+((
     $	C1(nu,nud,l)*dsqrt(dfloat(l+m)*dfloat(l-m)/
     $	dfloat(2*l-1))*psi000(nud,l-1,m)+			! NB this has changed since the calculations for makabe
     $	C2(nu,nud,l)*dsqrt(dfloat(l-m+1)*dfloat(l+m+1)/
     $	dfloat(2*l+3))*psi000(nud,l+1,m))/alf)
	
	tmp2=tmp2+((
     $	C1(nu,nud,l)*dsqrt(dfloat(l+m)*dfloat(l+m-1)/2.d0/
     $	dfloat(2*l-1))*psi000(nud,l-1,m-1)-
     $	C2(nu,nud,l)*dsqrt(dfloat(l-m+1)*dfloat(l-m+2)/2.d0/
     $	dfloat(2*l+3))*psi000(nud,l+1,m-1))/alf)
	endif
442	continue	

	 temporary=sqrootoftwo*zeta*
     $     (tmp1*alf-psi000(0,1,0)*psi000(nu,l,m))
	 rhs110(i)=real(temporary)
	 rhs110(i+noe)=dimag(temporary)
	 temporary=sqrootoftwo*zeta*
     $     (tmp2*alf+psi000(0,1,-1)*psi000(nu,l,m))
	 rhs11p1(i)=real(temporary)
	 rhs11p1(i+noe)=dimag(temporary)

	  i=i+1
450	continue	
470	continue	
500	continue	
		 
	call initial(aflag0,aflag1,iflag0,iflag1)
	call y12mbf(mmm,nz0,AA0,jcol0,ny,irow0,ny,ha0,nn,aflag0,
     +               iflag0,ifail0)
	if (ifail0.ne.0) go to 5152
	call y12mcf(mmm,nz0,AA0,jcol0,ny,irow0,ny,pivot0,rhs110,
     +       ha0,nn,aflag0,iflag0,ifail0)
	if (ifail0.ne.0) go to 5152
	call y12mdf(mmm,AA0,ny,rhs110,pivot0,jcol0,ha0,nn,iflag0,ifail0)
	iflag0(5)=3
 	call y12mdf(mmm,AA0,ny,rhs11p1,pivot0,jcol0,ha0,nn,iflag0,ifail0)

5152	continue
	
	i=1
	do 520 l=0,lmax
	 if (l.ge.mmax) then
	  	mmx=mmax
	 else
	 	mmx=l
	 endif	
	 do 520 m=-mmx,mmx
	   do 510 nu=0,numax
	   if ((nu.eq.0).and.(l.eq.0)) goto 510
       psi110(nu,l,m)=rhs110(i)+dcmplx(0.d0,1.d0)*rhs110(i+noe) 
       psi11p1(nu,l,m)=rhs11p1(i)+dcmplx(0.d0,1.d0)*rhs11p1(i+noe) 
	      i=i+1
510         continue
520      continue
	psi110(0,0,0)=0.d0
	psi11p1(0,0,0)=0.d0

	endif   ! of ndg.ge.0

c------ Calculating Transport Coefficients ---------------------------------

c...... Rection and Excitation rate ....................................................

	omegaonn_o=rr/dsqrt(muon2k/Tb)*sigma0	



c...... Average energy coefficients .....................................

	  tmpenergya=1.5d0*kb*Tb
	   tmpa1=tmpenergya*(real(psi000(0,0,0))-
     +	        dsqrt(2.d0/3.d0)*real(psi000(1,0,0)))
	   tmpa2=-tmpenergya*2.d0/dsqrt(3.d0)*real(psi11p1(1,0,0))
	   tmpa3=-tmpenergya*2.d0/dsqrt(3.d0)*dimag(psi11p1(1,0,0))
	   tmpa4=-tmpenergya*dsqrt(2.d0/3.d0)*real(psi110(1,0,0))
        en0=tmpa1/dabs(e)
	tmpenergy=dsqrt(muon2k/Tb)/sigma0/zeta/alf/sqrootoftwo
	gamx=tmpenergy*tmpa2
	gamy=tmpenergy*tmpa3
	gamz=-tmpenergy*tmpa4
	gam=dsqrt(gamx**2.d0+gamy**2.d0+gamz**2.d0)
c	theta=dacos(gamz/gam)*360.d0/2.d0/pi
c	phi=datan(gamy/gamx)*360.d0/2.d0/pi
	write(34,228) lmax,mmax,numax,gamx,gamy,gamz,gam,theta,phi
c	write(60,*) psideg,en0
c	write(61,*) psideg,gamx
c	write(62,*) psideg,gamy
c	write(63,*) psideg,gamz
c	write(64,*) psideg,gam
c	write(65,*) psideg,theta
c	write(66,*) psideg,phi

c.....  Drift velocity Components .................................

	 tmpdrifta=1/alf
	 tmpa1=tmpdrifta*real(psi000(0,1,1))*sqrootoftwo
	 tmpa2=tmpdrifta*dimag(psi000(0,1,1))*sqrootoftwo
	 tmpa3=tmpdrifta*psi000(0,1,0)
630     continue
	W_x=-tmpa1		! Flux drift velocities
	W_y=tmpa2
	W_z=tmpa3
	Wtot=dsqrt(W_x**2.d0+W_y**2.d0+W_z**2.d0)
	write(30,227) lmax,mmax,numax,en0,W_x,W_y,W_z,Wtot,omegaonn_o
	print*, 'Drifts (x,z), Energy= ',W_x,W_z,en0
	tmpa1=0.d0		! Bulk corrections
	tmpa2=0.d0
	tmpa3=0.d0
	do 631 nud=0,numax
	  tmpa1=tmpa1+Nbar(0,nud,0)*a(0,0,nud)*real(psi11p1(nud,0,0))
	  tmpa2=tmpa2-Nbar(0,nud,0)*a(0,0,nud)*dimag(psi11p1(nud,0,0))
	  tmpa3=tmpa3-Nbar(0,nud,0)*a(0,0,nud)*real(psi110(nud,0,0))	
631	continue  
	W_xb=W_x+tmpa1/alf/zeta/sqrootoftwo
	W_yb=W_y+tmpa2/alf/zeta/sqrootoftwo
	W_zb=W_z+tmpa3/alf/zeta/sqrootoftwo
	Wtotb=dsqrt(W_xb**2.d0+W_yb**2.d0+W_zb**2.d0)
	!write(32,227) lmax,mmax,numax,en0,W_xb,W_yb,W_zb,Wtotb,omegaonn_o
      write(32,227) lmax,numax,en0,W_z,W_zb,omegaonn_o
c	theta=dacos(W_z/Wtot)*360.d0/2.d0/pi
c	phi=datan(W_y/W_x)*360.d0/2.d0/pi
c	theta=0.d0
c	phi=0.d0

c	write(40,*) psideg,W_x
c	write(41,*) psideg,W_y
c	write(42,*) psideg,W_z
c	write(43,*) psideg,Wtot
c	write(44,*) psideg,theta
c	write(45,*) psideg,phi


c.....  Diffusion tensor components ..............................

	if (ndg.gt.0) then
	tmpdiffa=1/alf**2.d0
        tmpa1=tmpdiffa*(real(psi11p1(0,1,1))-
     $         real(psi11p1(0,1,-1)))
        tmpa2=tmpdiffa*(real(psi11p1(0,1,1))+
     $	       real(psi11p1(0,1,-1)))
        tmpa3=tmpdiffa*real(psi110(0,1,0))
        tmpa4=tmpdiffa*(dimag(psi11p1(0,1,1))-
     $	       dimag(psi11p1(0,1,-1)))
        tmpa5=tmpdiffa*real(psi110(0,1,1))
     $	       *sqrootoftwo
        tmpa6=tmpdiffa*(dimag(psi11p1(0,1,1))
     $	       +dimag(psi11p1(0,1,-1)))
        tmpa7=tmpdiffa*dimag(psi110(0,1,1))
     $     	*sqrootoftwo
        tmpa8=tmpdiffa*real(psi11p1(0,1,0))
     $          *sqrootoftwo
        tmpa9=tmpdiffa*dimag(psi11p1(0,1,0))
     $		*sqrootoftwo
	tmpdiffusion=dsqrt(muon2k/Tb)/sigma0/zeta/sqrootoftwo
	DXX=tmpdiffusion*tmpa1
	DYY=tmpdiffusion*tmpa2
	DZZ=tmpdiffusion*tmpa3
	DXY=tmpdiffusion*tmpa4
	DXZ=-tmpdiffusion*tmpa5
	DYX=-tmpdiffusion*tmpa6
	DYZ=+tmpdiffusion*tmpa7
	DZX=-tmpdiffusion*tmpa8
	DZY=-tmpdiffusion*tmpa9

	write(31,230) lmax,mmax,numax,DXX,DYY,DZZ,DXZ+DZX,DXY+DYX,DYZ+DZY
	write(33,230) lmax,mmax,numax,DXZ,DZX,DXY,DYX,DYZ,DZY
c	write(50,*)psideg,DXX
c	write(51,*)psideg,DYY
c	write(52,*)psideg,DZZ
c	write(53,*)psideg,DXY
c      write(54,*)psideg,DYX
c	write(55,*)psideg,DXZ
c        write(56,*)psideg,DZX
c	write(57,*)psideg,DYZ
c	write(58,*)psideg,DZY
c	write(80,*)psideg,DXY+DYX
c	write(81,*)psideg,DXZ+DZX
c	write(82,*)psideg,DYZ+DZY
c	write(83,*)psideg,DXX/DYY
c	write(84,*)psideg,DXX/DZZ
c	write(85,*)psideg,DYY/DZZ

	endif
c...... Calculation of the temperature tensor --------------------

	tmpa1=0.d0
	tmpa2=0.d0
	tmpa3=0.d0
	tmpa4=0.d0
	tmpa5=0.d0
	tmpa6=0.d0
	 tmpa=-sqrootoftwo*real(psi000(0,1,1))
	 tmpb=sqrootoftwo*dimag(psi000(0,1,1))
	 tmpc=real(psi000(0,1,0))
	     tmpa8=0.d0

     	 tmpa1=real(psi000(0,0,0))-dsqrt(2.d0/3.d0)*
     $		real(psi000(1,0,0))-dsqrt(1.d0/3.d0)*
     $  	real(psi000(0,2,0))+dsqrt(2.d0)*
     $		real(psi000(0,2,2))	

     	 tmpa2=real(psi000(0,0,0))-dsqrt(2.d0/3.d0)*
     $		real(psi000(1,0,0))-dsqrt(1.d0/3.d0)*
     $  	real(psi000(0,2,0))-dsqrt(2.d0)*
     $		real(psi000(0,2,2))	

     	 tmpa3=real(psi000(0,0,0))-dsqrt(2.d0/3.d0)*
     $          real(psi000(1,0,0))+dsqrt(4.d0/3.d0)*
     $  	real(psi000(0,2,0))

	 tmpa4=-sqrootoftwo*dimag(psi000(0,2,2))
	 tmpa5=-sqrootoftwo*real(psi000(0,2,1))
	 tmpa6=+sqrootoftwo*dimag(psi000(0,2,1))

	Txx=Tb*(tmpa1-tmpa*tmpa)	
	Tyy=Tb*(tmpa2-tmpb*tmpb)	
	Tzz=Tb*(tmpa3-tmpc*tmpc)	
	Txy=Tb*(tmpa4-tmpa*tmpb)	
	Txz=Tb*(tmpa5-tmpa*tmpc)	
	Tyz=Tb*(tmpa6-tmpb*tmpc)	
	Te=(Txx+Tyy+Tzz)/3.d0

	write(35,230) lmax,mmax,numax,Txx,Tyy,Tzz,Txy,Txz,Tyz,Te
c	write(70,*) psideg,Txx
c	write(71,*) psideg,Tyy
c	write(72,*) psideg,Tzz
c	write(73,*) psideg,Txy
c	write(74,*) psideg,Txz
c	write(75,*) psideg,Tyz
c	write(76,*) psideg,Txx/Tyy
c	write(77,*) psideg,Txx/Tzz
c	write(78,*) psideg,Tyy/Tzz
	
	!	DC 
c		Calculating excitation and power deposition rates


	do 9876 i=1,2
	 if(np(i).lt.0) np(i)=0
	 do 9876 j=-np(i),np(i)
	  tmp1=0.d0
	  tmp2=0.d0
	  do 9875 nu=0,numax   
		 tmp1=tmp1+vrate0(i,j,nu)*real(psi000(nu,0,0))
		 tmp2=tmp2+vrate1(i,j,nu)*real(psi000(nu,0,0))
c		 tmp1=tmp1+Nbar(0,nu,0)*vrate0(i,j,nu)*real(psi000(nu,0,0))
c		 tmp2=tmp2+Nbar(0,nu,0)*vrate1(i,j,nu)*real(psi000(nu,0,0))

9875	  continue
	  rate(i,j)=tmp1*sigma0/dsqrt(muon2k/Tb)
	  prate(i,j)=1.5d0*kb*Tb/e*(tmp1-dsqrt(2.d0/3.d0)*tmp2)*
     $				sigma0/dsqrt(muon2k/Tb)
9876	continue
	field_power=1.d-21*eonn*drift	


	
c	WRITE(17,*)(1.d0/nbar(0,j,0)*vrate0(1,14,j),j=0,4)

      write(36,539) lmax,numax,en0,W_z,W_zb,Dxx,Dzz,omegaonn_o,
     &  rate(1,0),rate(1,1),rate(1,2),rate(1,3),rate(1,4),Txx,Tzz,Gamz

!      write(36,227) lmax,numax,en0,W_z,W_zb,omegaonn_o
	!!!!!!!
	write(21,230) lmax,mmax,numax,rate(1,1),rate(1,2),rate(1,3),
     $ rate(1,14),rate(1,5)
	write(22,230) lmax,mmax,numax,rate(1,6),rate(1,7),rate(1,8),
     $ rate(1,9),rate(1,10)
	write(23,230) lmax,mmax,numax,rate(1,11),rate(1,12),rate(1,13),
     $ rate(1,14),rate(1,15)

C     $  rate(17),rate(12)+rate(13)+rate(14)+rate(15)+rate(16)
c	write(36,230) lmax,mmax,numax,rate(6),rate(8)+rate(9),rate(10),
c     $  rate(12),rate(7),rate(11)

	
c	write(21,230) lmax,mmax,numax,rate(1),rate(2),rate(3),rate(4),
c     &				rate(5) 
c	write(22,230) lmax,mmax,numax,rate(6),rate(7),rate(8),rate(9),
c     &				rate(10)
c	write(23,230) lmax,mmax,numax,rate(11),rate(12),rate(13),rate(14)
c     &				,rate(15)
c	write(24,230) lmax,mmax,numax,rate(16),rate(17),rate(18),rate(19)
c     &				,rate(20)
	
c	write(32,426) lmax,numax,rate(1,0,numax),rate(1,1,numax),
c     $ rate(1,2,numax),rate(1,3,numax),rate(1,4,numax),rate(1,5,numax)
c	write(33,426) lmax,numax,rate(1,6,numax),rate(1,7,numax),
c     $ rate(1,8,numax),rate(1,9,numax),rate(1,10,numax),rate(1,11,numax)
c	write(34,426) lmax,numax,rate(1,-1,numax),rate(1,-2,numax),
c     $ rate(1,-3,numax),rate(1,-4,numax),rate(1,-4,numax),
c     $  rate(1,-6,numax)
c	write(34,426) lmax,numax,prate(1,0,numax),prate(1,1,numax),
c    $ prate(1,-1,numax),prate(1,2,numax),prate(1,-2,numax),field_power
c	write(35,426) lmax,numax,prate(2,0,numax),prate(2,1,numax),
c    $ prate(2,-1,numax),prate(2,2,numax),prate(2,-2,numax),field_power


c	write(36,539) lmax,numax,en0,W_z,W_zb,Dxx,Dzz,Txx,Tzz,Gamz
539     format(1x,i1,1x,i2,14(1pd12.4))

998     continue                ! End of nu convergence test loop  
	write(30,*)
	write(31,*)
	write(32,*)
	write(33,*)
	write(34,*)
	write(35,*)
	write(36,*)
	write(21,*)
	write(22,*)
	write(23,*)
	write(24,*)

999     continue                ! End of l convergence test loop  

c------ Calculating Distribution Functions ---------------------------------
	
		dist0 = 0
		if (dist0.eq.0) then
	print*, 'in dist1'
        en=1.d-5
        pp=0                                    ! Energy counter
c        do 7000 k=0,6
	  k=0
          enptsi=enpts(k)
           do 7100 j=1,enptsi
            pp=pp+1
            
            energy(pp)=en
            c=dsqrt(2.d0*en*e/wi/u)     ! Speed
	    x=en*e/kb/Tb
            do 7200 l=0,1
	    ll=dfloat(l)
	     do 7105 m=0,l	
c	     m=0
	     dfm=dfloat(m)
                f0i=dcmplx(0.d0,0.d0)
	       
                do 7300 nu=0,nux

                  sv=dsqrt(dexp(faclog(nu+1)-gamlog(nu+l+2)))
     +           *sonine(nu,l,x)

                  f0i=f0i+sv*psi000(nu,l,m)
7300            continue
                svv=dsqrt(x**ll)*dexp(-x)/(kb*Tb/e)**1.5d0    
     $		     *(-1.d0)**dfm*dsqrt(dexp(faclog(l-m+1)
     $		     -faclog(l+m+1))*dfloat(2*l+1))	
     $		     	*dsqrt(2.d0)/pi**0.25d0

                f0(l,m,pp)=svv*f0i

7105	    continue	
7200       continue
c        write(92,*) energy(pp), dsqrt(energy(pp))*
c        write(92,*) energy(pp),
c     +				real(f0(0,0,pp))
cc        write(93,*) energy(pp), energy(j)*
c        write(93,*) energy(pp),
c     +                          real(f0(1,0,pp))
c	write(92,*) psideg,energy(pp),
c     +		dabs(real(f0(0,0,pp))) 
c	write(93,*) psideg,energy(pp),
c     +		dabs(real(f0(0,0,pp))) 
c	write(94,*) psideg,energy(pp),
c     +		dsqrt(energy(pp))*dabs(real(f0(0,0,pp))) 
c
	write(100,2270)  energy(pp),(dabs(real(f0(0,0,pp)))),
     +		(dabs(real(f0(1,0,pp)))),
     +  dsqrt(2.d0)*(dabs(real(f0(1,1,pp)))),
     +  dsqrt(2.d0)*(dabs(imag(f0(1,1,pp))))
     		
2270     format(5(1pd12.4))

			en=en+delen(k)

7100              continue

c7000            continue
        write(92,*)
        write(93,*)


	end if ! if dist0

c---  Plotting the velocity distribution function ----------------

	stop

	temp=dabs(dmod(psideg,vdfspacing))
	if ((temp.lt.1.d-9.and.psideg.ge.vdfspacing).or.
     +		(psideg.eq.0.d0)) then

        call assoclegendre(lx,0.d0,ALP) ! ie. theta=90.d0

	print*, 'counters=',jjj,vnn
	vnn=vnn+1	! Increment vdf counter

        noexp=enpts(1)
        noeyp=enpts(2)
	delenx=delen(1)
	deleny=delen(2)
	
	print*, 'phiplane=',phiplane
        do 7010 i=0,2*noexp
          enx=(delenx*(i-noexp))

cc	Calculation for phi plane
c	  if (enx.ge.0.d0) then
c	   phi=phiplane
c	  else
c	   phi=phiplane+pi 
c	  endif

          xen(i)=enx

          do 7020 j=0,2*noeyp
            eny=(deleny*(j-noeyp))
            yen(j)=eny
            en=dsqrt(enx**2.d0+eny**2.d0)
	    x=e*en/kb/Tb
	
c	Calculation in phi slice

c	   costheta=eny/en
c	   if (en.eq.0.d0) costheta=1.d0
c	    write(106,*) costheta
c           call assoclegendre(lx,costheta,ALP)

c	Calculation in theta=0.d0 slice

	   phi=datan(eny/enx)
	   if (enx.lt.0.d0) phi=phi+pi
	   if (enx.eq.0.0d0.and.eny.eq.0.d0) phi=0.d0

            disti=0.d0
            do 7030 l=0,lx
	    ll=dfloat(l)
	     do 7035 m=0,l	
	     dfm=dfloat(m) 
	 	f0i=dcmplx(0.d0,0.d0)
                do 7040 nu=0,nux

                  sv=dsqrt(dexp(faclog(nu+1))/dexp(gamlog(nu+l+2)))
     +           *sonine(nu,l,x)

                  f0i=f0i+sv*psi000(nu,l,m)
7040            continue

	        if (ll.eq.0.and.x.eq.0.d0) then
                svv=1.d0/(kb*Tb/e)**1.5d0    
     $		     *dsqrt(dexp(faclog(l-m+1))
     $		     /dexp(faclog(l+m+1))*dfloat(2*l+1))	
     $		     	*dsqrt(2.d0)/pi**0.25d0
	       else
                svv=dsqrt(x**ll)*dexp(-x)/(kb*Tb/e)**1.5d0    
     $		     *dsqrt(dexp(faclog(l-m+1))
     $		     /dexp(faclog(l+m+1))*dfloat(2*l+1))	
     $		     	*dsqrt(2.d0)/pi**0.25d0
		endif

		disti=disti+svv*(real(f0i)*dcos(dfm*phi)
     $			-dimag(f0i)*dsin(dfm*phi))
     $			*(-1.d0)**dfm*(2.d0-delta(m,0))*ALP(m,l)

7035	     continue
7030        continue
            distfn(vnn,i,j)=disti
 
7020      continue
7010    continue

          endif


	write(94,*) 

9999	continue		! Angle loop

	write(30,*)
	write(31,*)
	write(32,*)
	write(33,*)
	write(34,*)
	write(35,*)
	write(36,*)
	write(21,*)
	write(22,*)
	write(23,*)
	write(24,*)

	if (dist.eq.0) then
	write(91,*) vnn	
        write(91,*) 2*noexp+1
        write(91,*) 2*noeyp+1
	write(91,*) noexp*delenx
	write(91,*) noeyp*deleny
	write(91,*) delenx
	write(91,*) deleny
        do 561 i=0,2*noexp
          write(91,*) xen(i)
561     continue
        do 562 i=0,2*noeyp
          write(91,*) yen(i)
562     continue
	do 564 jjj=1,vnn
        do 563 i=0,2*noexp
         do 563 j=0,2*noeyp
          write(91,*) distfn(jjj,i,j)
563     continue
564	continue
	endif

227     format(1x,i1,1x,i1,1x,i2,6(1pd12.4))
228     format(1x,i1,1x,i1,1x,i2,4(1pd12.4),1x,0pf6.1,1x,f6.1)
229     format(13x,2(1pd15.7))
230     format(1x,i1,1x,i1,1x,i2,7(1pd12.4))
231     format(1x,i1,1x,i2,4(1pd12.4),1x,0pf6.1,1x,f6.1)
340     format(1x,f6.1,4(1pd13.5))
341	format(1x,f6.1,1pd13.5)
342     format(1x,i2,1x,i2,6(1pd14.6))


	end


	FUNCTION DELTA(I,J)

	INTEGER I,J
	DOUBLE PRECISION DELTA

	IF (I.EQ.J) THEN
		DELTA=1D0
	ELSE
		DELTA=0D0
	ENDIF
	RETURN
	END
	SUBROUTINE LEGP(L,X,P)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION P(10)
        P(1)=1.D0
        IF(L.EQ.1)GO TO 10
        P(2)=X
        IF(L.EQ.2)GO TO 10
        DO 5 N=3,L
        AN=DFLOAT(N-1)
        AN1=AN-1.D0
        AN2=2.D0*AN-1.D0
 5      P(N)=(AN2*X*P(N-1)-AN1*P(N-2))/AN
 10     RETURN
        END

	SUBROUTINE assoclegendre(lmax,x,p)

	integer lx
	parameter (lx=3)

	integer
     $  l,m,llimit
	double precision
     $  p(0:lx,0:lx),x,tmp,tmp1

c	Define low order p(m,l)

	p(0,0)=1.d0
	p(0,1)=x
	p(1,1)=dsqrt(1.d0-x**2.d0)
	p(0,2)=0.5d0*(3.d0*x**2.d0-1)
	p(1,2)=3.d0*x*p(1,1)
	p(2,2)=3.d0*p(1,1)**2.d0
	p(0,3)=0.5d0*x*(5.d0*x**2.d0-3.d0)
	p(1,3)=1.5d0*(5.d0*x**2.d0-1.d0)*p(1,1)
	p(2,3)=5.d0*x*p(2,2)
	p(3,3)=5.d0*p(1,1)*p(2,2)

c	write(104,*) x
c	do 5 l=0,3
c	  do 5 m=0,l
c	   write(104,*) l,m,p(m,l)
c	   p(m,l)=0.d0
c5	continue
c	p(0,0)=1.d0

c	if (lmax.gt.3) then

	tmp=dsqrt(1.d0-x**2.d0)
	do 10 l=1,lmax
	  tmp1=dfloat(2*l-1)*p(l-1,l-1)
	  p(l,l)=tmp*tmp1
	  p(l-1,l)=x*tmp1
10	continue

	do 20 m=0,lmax-2
c	  llimit=min(m+2,lmax)
	  llimit=m+2
	  do 30 l=llimit,lmax
	    p(m,l)=(x*dfloat(2*l-1)*p(m,l-1)-dfloat(l+m-1)*
     $			p(m,l-2))/dfloat(l-m)
30	  continue	
20	continue

c	endif
	
c	write(101,*) x
c	do 40 l=0,lmax
c	  do 40 m=0,l
c	   write(101,*) l,m,p(m,l)
c40	continue	

	return
	end

	FUNCTION sonine(n,l,x)

c	l = l+1/2 but we only pass l through since the arg of gamlog are
c	integers.  It is compensated for.

	INTEGER i,n,l
c	DOUBLE PRECISION x,m,sonine,faclog(200),gamlog(200),
c     +       vx(10,100,100),as(100,100,10)
 	DOUBLE PRECISION x,m,sonine,faclog(150),gamlog(150),
     +       vx(10,100,100),as(100,100,9)

c	DOUBLE PRECISION x,m,sonine,faclog(250),gamlog(250),
c     +       vx(10,200,200),as(200,200,9)
	COMMON/WEALTH/FACLOG,GAMLOG,VX,AS
	
c	(Gamlog(i)=dlog(Gamma((2*i-1)/2)) )
	       sonine=0.d0
		do 15 i=0,n
		  if ((x.eq.0).and.(i.eq.0)) then
		   sonine= sonine+
     + 		 dexp(gamlog(n+l+2))/dexp(gamlog(i+l+2))
     +  	 /dexp(faclog(i+1))/dexp(faclog(n-i+1))
	          else
	   	  sonine= sonine+
     +  		dexp(gamlog(n+l+2))/dexp(gamlog(i+l+2))
     +  		/dexp(faclog(i+1))/dexp(faclog(n-i+1))
     +        		*(-1.d0*x)**i
		  endif
15		continue
	return
	end
